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! ABSTRACT 

The X-ray binary LS I +61°303, consisting of a main sequence Be star and a compact 
object has been detected in the TeV range with MAGIC and VERITAS, and showed a 



clear intensity modulation as a function of the orbital phase. We describe a gamma-ray 
HjH attenuation model and apply it to this system. Our first result is that interaction of 

high energy photons with the background radiation produced by the main sequence 
Q H | star alone does not account for the observed modulation. We then include interactions 

between very high energy radiation and matter, and are able to constrain fundamental 
parameters of the system such as the mass of the compact object and the density of 
circumstellar matter around the Be star. In our analysis of the TeV data, we find that 
the compact object has mass M2 > 2.5M & at the 99% confidence level, implying it is 
most likely a black hole. However, we find a column density which conflicts with results 
from X-ray observations, suggesting that attenuation may not play an important role 
ON ' in the modulation. 

in 

O ■ 1. Introduction 



In the past few years, several high mass X-ray binaries have been detected as gamma ray emit 



ters (jAharonian et alj|200d : lAlbert J., et al.ll2006l ; lAcciari et al.ll2008l ). causing an intensification of 



observational and theoretical interest. High energy emitting binary systems consisting of a main 
sequence star and a compact object are the only known very high energy (VHE) galactic variable 
sources, and their short periods of days or weeks make them even more interesting observational 
targets. With increasing spectral coverage and statistics, the nature of photon emission and ab- 
sorption mechanisms is becoming increasingly constrained. Here we are concerned with high energy 
(TeV) photons emitted from the vicinity of the compact object and interacting with the background 
black body radiation and ejected material from the companion star. Even though these systems 
can be incredibly complex, a simple model of the absorption mechanisms and how they affect the 
system's light curve, can still shed light on many aspects such as the compact object mass and the 
orbital parameters. 
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One such example is the high energy emitti ng binary LS I + 61°3 03 ([Massi. et al.l 12004 ). It 
was first detected in the TeV range with MAGIC ([Albert J., et al.ll2006l ) and f urther observed wit h 
VERITAS at flux levels ranging between 5% and 20% of the Crab Nebula (|Acciari et alJbood ). 
This source has been observed throughout most o f the electrom agnetic spectrum starting with radio 
frequencies and extending to VHE gamma rays (jLeahvl 12004 ) . This broad spectral s tudy indicates 
that the system consists of a main sequence Be star of mass Mi = 12.5 ± 2.5 (Casares. el al 



20051 ) , surrounded by a circumstellar disk (iGrundstrom et al.1 120071 : IParedes et al.1 I20071 ). and a 
compact companion separated by tens o f solar radii at peria stron. The compact companion can 
be either a neutron star or a black hole ( Casare s. el al.ll2005l ). and its exact nature is still subject 



of investi gation and debate (IZdziars ki et al 



apastron (jAcciari et al. 



2008; 



Albert J., et al. 



201 ol) . The maximum VHE emission occurs close to 



20061 ) , suggesting that absorption plays an important 



role in the modulation. 



The outline of the paper is as follows: first we describe the model of attenuation due to pair 
production. Then we consider the particular case of LS I +61 303. We assume that the high energy 
radiation is emitted from the vicinity of the compact object and that its emission is isotropic and 
constant in time. The modulation due to photon-photon interactions is found to be insufficient to 
account for the VERITAS observations. For this reason, we include additional interactions of VHE 
photons with circumstellar material. This model permits to constrain the orbital parameters and 
the mass of the compact object as well as the density of ejected material from the companion star. 



2. Interaction with background radiation 
2.1. Radiative transfer equation 



The radiative transfer equation (| Chandr asekharl 1 1 96G ) for the intensity I(s,E), where s is the 
distance traveled by a photon of energy E from the emission point is 



dI(s,E) 
ds 



;i - cos n(s, e) a(E, e, £) I(s, E) + j(s, E) ; 



(1) 



Where n(s, e) is the spectral density of background photons of energy e emitted by the main 
sequence star, a(E, e, £) is the cross section^ for the interaction between photons colliding at angle 
£, and j(E, s) is a source term. 



Note that the term (1 — cos£(s')) corresponds to the relative velocity between the incident and target photons 
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2.2. Neglecting the source term 

The source term corresponds to secondary gamma-rays in the electromagnetic cascade. We 
can estimate the effect of the source term in the context of a cascade toy model, where the energy 
of these secondary gamma-rays is degraded by typically a factor of 4 (after the primary photon first 
produces an e + e~ pair). Consequently, the increase in the intensity I(s, E) of gamma rays of energy 
E is increased by 2tI(s,4E) + C(r 2 ), where r is the probability of interaction (or optical depth as 
defined in section [2T3|) . Assuming a power spectrum (section^]) we note that I(s, 4E) = I(s, E)4~^ . 
The source term is then approximately given by 



j(s, £)«2x 4-Tfi(s, e)a(E, e, £) I(s, E). 



(2) 



In this order of magnitude approximation, we have assumed that the cross section is the 
same for the ener gies E and AE. In the TeV range, the cross section decreases as ~ (ei?) _1 IneE 
(j Aharonianl 12004! ) . therefore eq. [2] constitutes an upper bound. Now taking 7 ~ 2, we note that 
the source term is smaller than the attenuation term in eq. [2 by a factor of < 10. Consequently, 
we ne glect the source term (For a detailed numerical calculation, see, e.g. ISierpowska fc Bednarek 

hood )). 



2.3. Solution of the Radiative transfer equation 

With the source term neglected, the solution to the radiative transfer equation is 



I(s,E)=I(s ,E) exp 



(1 - cos £(/)) n(s', e')a(E, e', s')ds'de' } . 



(3) 



Here so is the emission point in the vicinity of the compact object (see figured]), and e corre- 
sponds to the threshold energy for pair production, 



2mgC 4 



E(l -cosf(s))' 



(4) 



The dependence of the scattering angle £ in eq. [3] has been changed to a dependence on the 
path s. The problem then re duces to calculating the inte gral in the exponential of eq. [31 also known 
as the optical depth t(s,E) (IRibicki &: Lightmanlll979l ). In our calcu lation, we con sider the main 
sequence star as a point source, and in view of the results obtained bv iDubusI (120061 ). including the 
angular extension does not change our results significantly. 



The distribution of background black body photons can be taken as 
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n(r,e) = rc Q (e)-f , 



(5) 



where ro and no are the radius of the Be star and the density of background photons at this 
radius, i.e. 



27re 2 de /r \2 1 
nir ' £j ~ c 3 /i 3 e^-l' 

Here, the photon density has already been integrated over the solid angle. 



(6) 



3. The case of LS 1+61 303 
3.1. Attenuation 

There is debate as to what is the mechanism responsible for high energy emission. However, 
the aim of this paper is not to model the gamma ray emission but rather to investigate the effects 
of attenuation. This allows to derive a few characteristics of the main sequence star environment 
and compact object orbit. 



Grundstrom et al. (j2007l ) reported a temperature of T « 3 x 10 4 K and radius of R « 6.7i? G 



■0 

for the Be star. The black body distribution peaks at a few eV, and the threshold energy for pair 
production with TeV incident photons is of the order of 1 eV, so that most of the background 
photons may contribute to the attenuation, provided the scattering angle is favorable. The back- 
ground photon density is found to be of the order of n 7 



~ 10 12 cm 3 at a the radius of the star. 



The ci rcumstellar disk has been observed by Waters et al. ([Waters et al.lll9 88) and by Paredes 
et al. (jParedes et al.l 120071 ). who estimate the disc ion density to be n e ~ 10 13 cm -3 at one stellar 
radius. The cross section for pair production is of the order of <r 77 « O.lor at the threshold energy. 
The cross section for interaction w i th hydrogen has a constant value of <t 7 # ~ 2 x 10~ 2 <7t above 



a few hundred MeV (Heftier 



1954 



Aharonianl 120041 ) . With these cross sections, a first estimate 



suggests that both interactions may result in comparable degrees of attenuation. However, there is 
a strong angular dependence in the 77 interaction, the extreme case being when the both photons 
are emitted in the same direction, a configuration in which there is no VHE attenuation. Also the 
threshold energy is much higher when the incident and target photons are nearly parallel, so fewer 
background photons contribute to attenuation. Consequently, 77 attenuation may not have strong 
modulation as a function of the orbital phase when compared with the modulation produced by 
interactions with the circumstellar material. 
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3.2. Orbital parameters of LSI + 61°303 



The orbital parameters of LS I + 61°303 , as sketched in figur e [Q are still subject of research 



(jAragona et al 



2009 



Grundstrom et al.l 120071 ; ICasares. el al.ll2005l ). The quantities of interest are: 



the period P, the angle between the major axis of the ellipse and the line of sight ip, the projected 
semi-major axis (ai sinz), corresponding to the ellipse of the Be stai@, the eccentricity rj, the phase 
at periastron 0o, and the mass function /(Mi,M2), which depends on the period and the radial 
velocity and relates the masses of both objects and the inclination angle i. The most recent 
orbital solution has been obtained by lAragona et al.l (|2009l ). where P = 26.4960 d, ip = 40.5 ±5.7°, 
aisini = 8.64 ± 0.52 R & , rj = 0.54 ± 0.03, <fo = 0.275 and f(M 1 ,M 2 ) = 0.0124 ± 0.0022 M . 
Assuming the attenuation is responsible for the periodic modulation of the TeV emission, our 
results constrain i and thereby the mass of the compact object. 




Fig. 1. — Sketch of the orbital parameters of LS I + 61°303 for the case of i = 90°. The arrow 
points to the observer. Also shown are the values of the phase at periastron, apastron, and superior 
and inferior conjuction. 



4. The integrated flux 



Following observations of LS7 + 61°303 from 09/2006 to 02/2008, the VERITAS collaboration 



reported power law spectrum = <I>o (j§y) 7 ) with a spectral index of 7 = 2A±0.2 K t n t ± 0.2 
at energies above ~ 0.5 TeV, and between phases (ft = 0.6 and <f> = 0.8 (jAcciari et al.l 12008). 



Observations at lower energies made by Fermi between 08/2008 and 03/2009, indicate that the 



2 The projected semi- major axis of the ellipse described by the compact object is typically labeled as 02 sin i. 
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spectral index does not change significantly as a function of the orbital phase (jAbdo et al.l 12009 ) 



Therefore, we _ a constant intrinsic^ spectrin as a taction of the phasT^^" . 
The integrated flux is then 

^ = 12 iSL 1 ^ ^ dE = F " £ (I) " 7 1 *><^ m 

where -Eo depends on the detection threshold energy of the detector, and Fq is a normalization 
factor that is taken as a free parameter. 



4.1. Light curve assuming only 77 interactions 

Figure [2] shows the attenuation as a function of the orbital phase for several different energies 
for the case of the compact object having the canonical neutron star mass (i ~ 64° or M ~ 1.5M & ). 
In figure [2] we essentially reproduce one of the results obtained by Dubus, except that the orbital 



parameters used are the newer set obtained bv lAragona et al.l (|2009l ). When only interactions with 



the background black body photons are taken into account, and the orbital plane is closer to being 
seen edge-on, the optical depth approaches a minimum when the compact object is close to the 
main sequence star. This is especially the case for very high inclination angles, corresponding to 
the mass of the compact object being small, and close to the Chandrasekhar mass. This behavior 
can be understood from the angular dependence of the threshold energy in addition to the relative 
velocity of the incident and target photons approaching a minimum. Also, at high energies, the 
cross section for pair creation decreases as the inverse square of the center of mass energy, decreas- 
ing the optical depth even more. That is, even though the total density of background photons 
increases (as 1/r 2 ) when the compact object approaches the Be star, a combination of the previ- 
ously mentioned factors dominates as can seen in figure[2l 



Figure [3] shows the normalized integrated flux assuming different incli nation angles and cor - 
responding compact object masses. The VERITAS data shown in figure [31 ( Weinstein et al.l 12008) 
were binned to show a single light-curve as opposed to monthly data. If we assume that the emission 
comes from the vicinity of the compact object, and is isotropic, and constant as a function of the 
orbital phase, then these results lead us to conclude that there must be an additional attenuation 
mechanism at play. 



3 By intrinsic we mean non attenuated by pair production. 
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Fig. 2. — Attenuation e~ T77 as a function of the orbital phase for different primary photon energies 
(77 interactions only). A mass of 1.5Mq, corresponding to i = 64°, was assumed for the compact 
object. 



4.2. Light curve including 77 and interactions 



The detailed structure of the circumstell ar material su r round i ng a Be star in t h e pre sence 
of a compact com panion has been studied by I Waters et al.l ()1988l ) , iMarti fe Paredes! (|1995l ) and 
Reig et al.l (|2000l ) among others. It is thought to consist of a main equatorial disk-like component 
and a polar wind. Typically, the parameters that describe the decretion disk include: the mass loss 
rate, the wind termination velocity, the half opening angle of the disk, and radius of the disk. The 
polar wind is radiatively driven and may have a corresponding m ass loss rate comparable to the 
equatorial wind (within one order of magnitude (IWaters et al.lll988l )). When comparing the quality 
of the data shown in figure El and the complexity of the models that describe the circumstellar 
material, it appears that, when only TeV data is used, only an order of magnitude estimate of the 
material density and its extension in the system can be achieved. With this in mind, we rather 
assume a simple isotropic distribution of material that decreases as a power q of the distance from 
the Be star (n = n#(ro/r) 9 ). We start by setting q = 2 and then consider different values for com- 
parison. Parameters found from existing models are taken into consideration for our approximation. 



For the case of a constant cross section and a 1/r 2 distribution of hydrogen, the optical depth 
can actually be found analytically (see appendix A). The resulting light curves are shown in figure 
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DC 




Fig. 3. — Normalized light-curve for 77 interactions only. Each curve corresponds to a different 
mass of the compact object and inclination angle. 



ID for different values of the inclination angle. For low inclination angles (high mass), VHE photons 
emitted from apastron go through less circumstellar hydrogen than those emitted from periastron. 
This results in minimum attenuation near apastron, as can be seen in figured! From this figure it 
is clear that the emission peak corresponding to a canonical 1.5M© neutron star is only marginally 
supported by observations. 



We take the inclination angle, characteristic density, and normalization factor as free parame- 
ters. We find the inclination angle to be i < 28° (M2 > 3M Q ) at the 89% confidence level (CL) in 
the context of this model, and % < 34° {M% > 2.5Mq) at the 99% CL. These limits are not in good 
agreement with the neutron star scenario generally favored for the broad-band spectrum it im- 
plied However, our results are still consistent with other observational constraints (10° < i < 60°) 
(jCasares. el al.ll2005l ) obtained from optical spectroscopy. As for the circumstellar material, if we 
assu me the characteristic extension to be tq ~ 100-Rq, consistent with more sophisticated mod- 
els (jSierpowska-Bartosik &: Torres! |2009| ) . then the density of hydrogen in the disk is found to be 



2.0 x 10 13 cm~ 3 < n H < 1.9 



10 15 cm 3 at the 



CL and n# 



x 10 14 cm 3 at the 



4 See Zdziarski et al. ( Zdziarski et al. 



2010) for more details 
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Fig. 4. — Light curves for different inclination angles i and hydrogen densities nn when an isotropic 
distribution of hydrogen is included. The mass of the compact object is 4.0M Q for the two curves 
whose peak is at phase ~ 0.7, and 1.5M© for the curve whose peak is closer to phase ~ 0.4. The 
concentration of hydrogen at r « 100i? Q for each curve is labeled in the top right-hand corner of 
the figure. In the two right-most curves, the intrinsic luminosity (normalization factor) is chosen 
so as to maximize agreement with data, and it increases with the hydrogen density. 
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CL. 



By integrating the volume density along the line of sight to the compact object at apastron, 
we find a column density of 1.9 x 10 26 cm -2 < Nh < 1. 8 x 10 28 cm~ 2 at the 99% CL, which is much 
higher than results f ound elsewhere in the literature (| Waters et al.l Il988l : iMarti fc Paredesl Il995l : 



Esposito et al. 1120071 ) . In part icular, when we use th e column density found by X-ray observations 



Nh = (5.7 ± 0.3) x 10 21 cm~ 2 lEsposito et al. I (|2007n . we find a reduced x 2 of 3.06 (11 degrees of 
freedom), corresponding to a \ 2 probability P(x 2 > 3.06) = 0.04%. A rough estimate suggests 
that by including ~ 10% of helium, the column density would be reduced by a factor of ~ 2, which 
is not sufficient to achieve compatibility with X-ray results. 



Density profiles in Be stars typically have radial dependences of l/r q , where 2.3 < q < 3.3 
(jLamers &: Waters! 1 198 71 ). depending on the opening angle of the disk. Therefore we expect our 
constraint on the density to constitute a lower bound_|. We perform our calculation with q = 3 and 
note that our results do not change considerably. 



The hydrogen density also corresponds to a mass loss rate of Mi ~ 10 
where O is the solid angle. Typically accepted values for the m a ss los s rat e are in the range o f 
~ lO^Moyr -1 to lO _8 M yr _1 , as have been reported by ISnowl ()198ll ) and Iwaters et al] (|1988l ) 
among others. A first glance at our result for the mass loss implies that it doe s not agree with 
the observations, i.e. setting f2 = An and V w i n d ~ 100km s -1 ( Waters et al.lll988l ). However, if we 
relax the assumption of an isotropic distribution of hydrogen, our result implies that small solid 
angles are favored as well as small velocities for the stellar wind. Small solid angles are consistent 
with the thin disk scenario that is most commonly accepted. Small velocities of the order of a few 



km s 1 are however not consistent with what is found elsewhere in the literature, e.g. ( Waters et al 



19881 ). and the wind indeed has higher velocities, this would imply that the system may have been 
observed while in a state of high mass loss rate. 



5. Discussion 

Since, in the TeV range, the interaction with matter is approximately independent of the 
energy, and since, as figure [3] shows, 77 interactions are insufficient to account for the orbital 
modulation, then the intrinsic non-attenuated differential spectrum is essentially the same as the 
observed spectrum (a power law of spectral index —2.4). However, the intrinsic TeV luminos- 
ity is several orders of magnitude higher than the measured luminosity. Taking the distance to 



This is assuming that the disk and orbit lie in the same plane. 
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the source to be approximately 1.8 kpc ([Steele et al.lll998l ). we find the intrinsic luminosity to be 
L ~ 5 x 10 37 ergs _1 when the hydrogen densi ty is of th e orde r of ~ 5 x 10 13 cm -3 . This intrinsic 
luminosity is comparable to that suggested bv lBottcher I (|2007l ) for LS 5039, the only other known 
TeV binary thought to contain a black hole. 



It is interesting to compare this intrinsic luminosity to the Eddington luminositjo L^dd ~ 
1.3 x 1O 39 (M2/M0)ergs _1 , which is comparable to L, and implies that radiation may be beamed 
in our direction. It is also interesting to calculate the accretion rate that would be needed in or- 
der to obtain the intrinsic luminosity: By taking L ~ GM2M2/R, where R is of the order of the 
Schwarzschild radius (2GM2/C 2 ), we find M2 ~ 2 x 10 il^yr . This rate is comparable with the 
observed mass loss rate of ~ lO~ 8 M yr _1 . The fact that the accretion rate is comparable to the 
measured mass loss rate, suggests that the flow of matter can be quite complicated, e.g. an increase 
in the accretion rate would strip most of the circumstellar mass, leading to time variability. This 
may explain the fact that no VHE detections have been reported since 2008. 



Still assuming the intrinsic luminosity to be constant in time, we can estimate the amount of 
hydrogen needed to attenuate the source to below the detectability threshold. We find that the 
density must increase from ~ 5 x 10 13 cm~ 3 to ~ 5 x 10 14 cm~ 3 at the characteristic distance of 
100i?Q. This amount of hydrogen in turn leads to much higher mass loss rates than those observed, 
and it may also imply a stronger activity of the source. 



It is worth mentioning that the attenuation model is not the only possible way to account 
for the modulation. For example, there is also the possibility of the emission being anisotropic, 
and the modulation r esulting from a geometrical effect. This possibility is described in detail by 



Zdziarski et al.l (|2010l ). were a shocked pulsar wind with a large Lorentz factor is thought to be the 



cause of emission. 



6. Conclusions 



For the case of LS I + 61°303, we find that attenuation due to 77 interactions with the back- 
ground radiation does not account for the observed high energy flux modulation as a function of 
the orbital phase, namely a narrow peak near apastron. This effect leads us to investigate some 
properties of the ejected material from the Be star, and the inclination angle of the orbit. We 
find the angle of the orbit to be i < 34° (M > 2.5M & ) at the 99% confidence level, suggesting 
that the compact object is a black hole rather than a neutron star. We also find the density of 
hydrogen in the disk to be 2 x 10 13 cm~ 3 < n H < 2 x 10 15 cm~ 3 at the 99% CL (at lOOi? ), which 



At the energies considered here, the cross section for inverse Compton is ~ O.lor 
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accounts for most of the observed gamma ray absorption. If the compact object is indeed a black 
hole as our analysis suggests, then the gamma ray emission is likely to be powered by accretion 
(jBosch- Ramon et al.l 120061 ; IZdziarski et al.ll2010n . Also, a black hole scenario might be even more 
complicated due to the possibility of VHE emission originating from termination of jets, therefore 
we cannot exclude the possibility of the modulation being due to geometrical effects. Current VHE 
data does not allow to constrain the system much more than what we ha ye already done, an d the 
fact that VHE detections hav e not been reported since the VERITAS (jAcciari et al.1 120081 ) and 
MAGIC ([Albert J., et al.ll2006l ) detections where made, makes the problem even more puzzling. A 
possible explanation might originate from a complex matter flow. This is suggested by the fact 
that the accretion rate needed to explain an intrinsic non-attenuated luminosity, is comparable to 
the measured mass loss rate of the Be star. 



An inconsistency arises when comparing our results with those derived from X-ray observations. 
We find the column density to be 1.9 x 10 26 cm~ 2 < N H < 1.8 x 10 28 cm~ 2 (99% CL), which is only 
compatible with X-ray results at the 0.04% confidence level. Such an incompatibility may imply that 
pair production in the stellar wind is not the cause of the modulation. Consequentl y, our estimates 



on the mass and column density may not be valid. An alternative explanation by IZdziarski et al 
(|2010l ) suggests that the modulation is due to a geometrical effect. Here a shocked pulsar wind is 
thought to flow along a cone with a large Lorentz factor, producing beamed radiation which can 
be seen when the cone passes through the line of sight. 



Appendix A: Optical depth for constant cross section and 1/r 2 density distribution 

Using a 1/r 2 distribution of hydrogen, the cross section an accounting for interactions between 
VHE photons and hydrogen, and the system of coordinates shown in figure[5] (corresponding to an 
orbital plane seen edge on) , we can calculate the integral for the optical depth to be 



°° n H rla H 



x 2 + yf + zf 



dx 



'yf + zf 



(8) 



riHTQCTH 



tan 



'yf + zf 



(9) 



where Xj, y^, and Zi are functions of the orbital angle 6, and rg is the characteristic radius of 
the hydrogen disk. For the case of a circular orbit as seen edge on (figure E|) , we can easily see 
the limiting behavior of the intensity as a function of the orbital angle. That is, expanding around 
8 ~ reveals that the attenuation around this region behaves like a Gaussian. 
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For 9~0: I(n, 9) = I (9, n)e ^ ( + '. (10) 

Similarly, expanding around 9 ~ tt/2 reveals that the attenuation behaves like a decreasing 
exponential 

For 9^ tt/2: I(n, 9) = I (9, n)e r * (11) 

For a more complicated geometry of LS I + 61°303, it is now just a matter of inserting the 
appropriate expressions for Xi{9), yi{9) and Zi{9). 
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Fig. 5. — Coordinate system used for calculating the optical depth (eq. [9]). 



